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A METHOD OF ESTIMATING THE CARRIER FREQUENCY OF A PHASE- 
MODULATED SIGNAL 

5 

The present invention relates to a method of estimating the carrier frequency of a phase- 
modulated signal. 

10 Modulation is the process of superimposing the characteristics of a signal onto a carrier 
wave so that the information contained by the signal can be transmitted by the carrier 
wave. In a phase modulated signal, the relative phase of the carrier wave is varied to 
encode the information contained in the signal. In communications it is often necessary 
to measure the frequency of the carrier wave. For most purposes this can be achieved by 

15 examining the frequency profile of the transmitted signal in the Fourier domain, where 
the- highest peak of the profile represents the strongest frequency (i.e., the carrier) 
component. In certain situations, however, the accuracy of phase measurements required 
is far greater than that obtained using a Fourier transform, and therefore more accurate 
techniques are required. 

20 

The phase 0(f) of a signal can be represented graphically by a vector in the real- 
imaginary plane, as shown in Figure 1. The amplitude of the signal is given by the 
magnitude of the vector, and the phase is given by the arc tangent of the ratio of the 
signals in the Q (real) and I (imaginary) channels, i.e., tan" 1 ^. Unfortunately 

25 individual phase measurements represented in this way are limited to within a range of 
size of 2n, as illustrated by the following example. The phase of a sinusoidal signal 
x(r) = sin(cur) increases linearly with time, and the gradient of the plot of this function is 
proportional to the frequency, go. However, since direct phase measurements are limited' 
to within the range of - n to + n, the phase plot appears as a "sawtooth" shape over time: 

30 initially the phase increases linearly until it reaches the value of -5-7t, it then jumps 
discontinuously to -n before increasing linearly again. The process of identifying and 
removing these discontinuities to recover the true linear phase relationship is known as 
phase unwrapping. Removing a discontinuity in the phase plot involves placing points 
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in the correct phase range by adding ±2iz to all points after the discontinuity, depending 
on whether the discontinuity is a positive or a negative jump. 

An example of a phase wrapped signal and the respective phase unwrapped signal is 
shown in Figure 2a and Figure 2b, respectively. It is obvious from Figure 2a that when 
trying to unwrap the wrapped phase, the problem lies in detecting the discontinuities. 

For a simple unmodulated carrier wave of constant frequency, an accurate measure of 
frequency can be gained by unwrapping the phase to produce a phase. plot, and fitting a 
straight line to this phase plot. The gradient of this straight line is proportional to the 
frequency of the signal. This technique is used in US Patent No. 5,019,823 (Thorn 
EMI). Accurate frequency measurements can be obtained using this technique because 
all the data in the signal are being used. However, when the data are noisy it is difficult 
to unwrap the phase to a straight line because phase discontinuities can occur over 
several samples and it can be difficult to correctly identify these phase discontinuities. If 
just one point is incorrectly placed in the wrong phase range the accuracy of the result is 
ruined, as each following point is unwrapped using the preceding point. 

For complicated signals, such as phase modulated signals, phase unwrapping becomes 
unreliable because of the presence of phase transitions. During phase transitions the 
frequency and phase of the signal varies rapidly, so it is difficult to identify the phase 
discontinuities and shift the proceeding phase values into the correct phase range. Even 
if phase unwrapping is performed correctly, the fact that the phase during the transitions 
covers the whole range [-71. .7t] stops a straight line from being fitted to the data. The 
conventional way around this is to know where the phase transitions are beforehand and 
then remove them. This allows the data to be unwrapped and a straight line to be fitted, 
but fails when noise is present. Of course, the situations where little noise is present and 
the positions of all the transitions are known to begin with are few and far between! 

A method is therefore required that neither relies on prior knowledge of the previous 
sample's phase (and is therefore robust to noise), nor requires the positions (or the 
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number) of transitions to be known. The method should also be independent of the exact 
form of the signal during phase transitions for widest applicability. 

An aim of the present invention is to provide a method for estimating the frequency of 
5 the carrier wave of a phase-modulated signal. A further aim of the present invention is 
to provide a method for estimating the frequency of the carrier wave of a phase- 
modulated signal in the presence of noise. 

An embodiment of the invention will now be described, by way of example only, with 
reference to the accompanying Figures, in which: - 

10 Figure 1 is a graphical representation of the measured phase <p(t) of a complex signal; 

Figures 2a and 2b show a graph of a wrapped and unwrapped signal, respectively; 

Figure 3 shows a graph of the real part of the complex signal; 

Figure 4 shows a graph of the phase of the raw data; 

Figure 5 shows a graph of phase data wrapped from 0 to n; 

15 Figure 6 shows a graph of the phase of the signal after one iteration of the method; 

Figure 7 shows a graph of the phase of the signal after a number of iterations of the 
method; 

Figure 8 shows a graph of the final phase plot; 
Figure 9 shows a block diagram of the main stages of the method; 
20 Figure 10 show a block diagram of a number of method steps; and 
Figure 1 1 shows a block diagram of the remaining method steps. 



The estimation of the carrier frequency f e of a signal with, for example, two phase- 
modulation states with added noise will now be described. 
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Firstly, in step one of the method, the signal pulse-train must be extracted from the 
signal (110). There are two ways this is commonly achieved: 1) a pulse-train is assumed 
to be a valid signal when its power is greater than a pre-determined level for the system, 
or 2) the use of a correlation technique with a template to find the start of the pulse train. 
5 However, knowing the exact starting point of the signal is not important to the outcome 
of the method. 

The second step of the method requires that an initial estimate f e % of the carrier 

frequency f c of the signal is made (112). This can be found from the highest peak of a 

Fourier transformation of the signal, i.e., in the frequency domain. This is a common 
10 signal processing technique and will therefore not be discussed in detail. 

Step three of the method attempts to mix the signal down to 0 Hz by removing the initial 
estimate f c " of the carrier frequency (120), resulting in complex (i.e., IQ) data. This 

procedure is well documented in standard signal processing texts, but is essentially 
carried out as follows. The signal is multiplied in the time domain by the complex 
15 function e l0a , where Q) = 27tf c ' and/ c ' is the estimated carrier frequency. This 

procedure transforms the signal from real to complex. The resultant signal has a 
component at the difference between f c ' and f c and an unwanted component at / c ' + 

f c which results from the well-known relationship 

sin(cy,r)*sinr(G) 2 r) = sin([a>, - w 2 ]r) + sin([<D 1 + co 2 ]t) where, in this case, co l and 

20 co 2 correspond to 27if c and 27tf c \ respectively. The high frequency component of the 

resultant signal is then removed using a low-pass filter. These steps are summarised by 
the following equation: 

l{t) = lowpassfilter\ s(t) cos(27if c ' ) - 
Q(t) = lowpassfHter-s(t)sir\{27rf c ') 

where s(t) is the phase modulated signal. I(t) is shown in Figure 3. From this Figure it 
25 can be seen that the initial estimate f e % of the carrier frequency was not accurate. If it 
had been accurate, only the phase transitions would remain. 



The iterative part of the method starts at step four. The first stage of this step is to 

calculate the phase 0(0 of the processed data (122), where <p(t) = tan '^ /( /^( r )] " The 

result of applying this stage of step four to the signal is shown in Figure 4. Two features 
are clearly seen from this Figure. Firstly, the phase increases with time until it reaches 
5 the value n, at which point the phase jumps discontinuously to -n before continuing its 
upward trend. Secondly, four phase transitions can be seen which interrupt this pattern. . 
It is these phase transitions that prevent conventional unwrapping techniques from being 
used. It should also be noted that the gradient of the phase versus time plot may be 
positive or negative, depending on whether the estimate f c ' of the carrier frequency is 

10 above or below the actual carrier frequency f c . 

The next stage of step four is to calculate the mean phase (124) and then to wrap the 
data into a phase range (126) equal to the size of the phase transition, making sure that 
the cyclic mean of the data is centred halfway through the range. As this example uses a 
15 signal with two phase states, we wrap the data into a range from 0 to k, and shift the 
cyclic mean to n 12. The wrapping equations used to perform this operation can be 
written as follows: 



20 



25 



shifted S 



<P(t) wrapped = 0(' W«, " ?r x floor 



n 



(1) 



where (p(t) is the phase, n is an integer large enough such that (j>(t) + niz is always 
positive, (0)is the cyclic mean of <j>(t) , and the function floor(x) takes the value of the 
largest integer smaller than x. The factor of 7t/2 is introduced because it is equal to half 
the size of the phase transition. 



When adding multiples of the phase jump (i.e., nil), the multiple must be large enough 
to make the phase positive. It should also be large enough to offset any large negative 
value from the cyclic mean and to offset the fact that (depending on the size of the phase 
transitions) the phase may be several transition multiples below zero. There is no 
30 maximum size for this value although if it is too large, computation accuracy may be 
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compromised. The cyclic mean can be in the range [0...2n] or [-7t...7r] since it will be 
wrapped down to the size of a phase transition anyway. Subtracting the cyclic mean and 
then adding half the phase range in Equation (1) shifts the mean phase of the wrapped 
data to the centre of the phase wrapping range. Phase unwrapping is achieved in this 
process by subtracting the integer number of multiples of the phase wrapping size from 
the shifted phase result. This ensures that the final range of the wrapped data is between 
0 and the size of one transition. The effect of wrapping the phase in this way is that 
when the gradient of the phase versus time plot is small enough, a phase transition in the 
original signal causes a jump in the wrapped phase that will start and end at the same 
phase value. 

The cyclic mean is defined as the mean phase of a group of complex data (such as IQ 
points in a communications signal) whose phase range wraps around from 0, to <p 2 , e.g. 
from -7r to 7r. The wraparound causes a simple mean of the phase angles to be invalid, 
which is the reason for using the cyclic mean. The cyclic mean can be visualised by 
plotting the phase vectors of the complex data on the Argand plane and then calculating 
the vector mean. The phase angle of the vector mean is the cyclic mean of the set of data 
points, and is calculated as follows: 



where N is the number of data samples, and 0. is the phase of the z'th sample. 

The reason for the use of the cyclic mean in this method is to move the maximum 
density of data points to the centre of the vertical range of the wrapped signal. Then, 
when the gradient of the phase versus time plots has been completely removed, the 
resulting horizontal trend line is in the centre of the range. If this were not done, it could 
result in the first half of the data disappearing off the top of the range, and the second 
half of the data reappearing at the bottom during an iteration of the algorithm. It would 
then be impossible to fit a line longer than one of the two segments of data, and 
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therefore the gradient estimation would not get any more accurate. Also, because the 
maximum density of data points is in the centre of the range, the outlying points during 
the phase changes are distributed evenly above and below the trend line. This minimises 
the effect of phase transitions on the gradient of the line of best fit to the data points. 

5 

As already stated, if data is centred on the cyclic mean, then the maximum data density 
is in the centre of the range of data points. This explanation has little intuitive meaning. 
The position of the mean is simple to understand if we consider a physical example. 
Imagine the graph shown in Figure 5 rolled into a tube, with the phase axis forming the 

10 circumference of the tube. Then, if each data point on the graph is represented by a 
small weight stuck to the inside of the tube, and the tube set down on a flat surface and 
allowed to roll freely, when it comes to rest the cyclic mean will be the line directly 
under the tube where it touches the surface. This is the point of greatest data density. If 
we further imagine the phase axis of the graph to be printed on a transparent sleeve 

15 covering the tube, we can rotate the sleeve so that the centre of the scale (x 12 here) is 
aligned with the cyclic mean at the point where the tube touches the flat surface. If the 
phase axis of the sleeve is now fixed relative to the data and the tube is unwrapped to a 
flat sheet once more, we will have a graph where the data is centred on the cyclic mean. 
The result of wrapping the data shown in Figure 4 using the method of the invention is 

20 shown in Figure 5. 

In step five of the method, short straight lines are fitted to the phase plot (i.e., to sections 
of H0 W rapped ) usin S a least-squares fit or similar technique (128). In the first iteration of 
the algorithm, it should be ensured that some of the lines are fitted to the non-transition 
25 sections of the signal. For this reason, the lengths of the straight lines are chosen to be 
substantially shorter than the expected distance between phase transitions. 

In the sixth step a x~ measure, for example, is used to examine how well the data fit a 
straight line model. The lines that fit well will be the ones that do not fall across phase 
30 transitions (which are not in general linear) or wrapping discontinuities. An average of 
the gradient (g) of these lines is taken (130). This average gradient is proportional to 
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the difference between the actual carrier frequency f c and the . estimated carrier 
frequency f c \ and is hereinafter referred to as the frequency error //. Another way of 
looking at this technique is to imagine we are trying to remove the average gradient of 
the phase versus time data plot so that it becomes a horizontal line (i.e. it corresponds to 
OHz). 

In the next step of the method, an improved estimation of the carrier frequency of the 
signal is given by adding the current approximation f c ' of the carrier frequency to the 
frequency error f e calculated in the previous step. The improved estimate of the carrier 
frequency must now be removed from the data. Two methods of achieving this are: 

• Mix down the original sample signal by the improved estimation of the carrier 
frequency using conventional signal processing methods, and reformulate the phase 
plot to the range 0 to n. This method has the advantage over other methods in that it 
is less susceptible to round-off error, although it tends to be computationally 
intensive. 

• Skew the data by the average gradient (140) using <p(t) = <p(t) - (g)t thereby trying 
to produce a horizontal line. Keep a running total of the skews that have been 
applied to the data. This method removes the need to calculate the sines and cosines 
required for conventional mixing down, and is therefore more suitable for real-time 
implementation as it is more efficient. Once we have reached the end of the 
algorithm, the sum of the average gradients measured at each iteration can be 
converted into a frequency which represents the quantity we are trying to measure, 

Using either. of these means of removing the current estimate f c ' of the carrier 
frequency results in the phase plot shown in Figure 6. 

The improved carrier frequency estimate f r ' is then calculated (142) using the 
following equation: f c '= f c '+j-^ . The line length / is then increased. 
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Method steps 4 to 7 are now repeated, this time fitting longer straight lines to the phase 
plot. The next estimate / c 'of the carrier frequency will be made from a larger number of 
data points than the initial estimate, and will therefore be more accurate. The length / of 
the lines fitted to the data is (usually) increased for each iteration of the algorithm until 
5 the line length extends the whole length of the signal. The increase in line length / after 
each iteration depends on the form of the signal and is likely to be application 
dependent. As a rule of thumb / can be doubled after each iteration. If / is increased too 
quickly between iterations then convergence on the optimum carrier frequency 
measurement cannot be guaranteed. This can be avoided by examining the proportion of 

10 "well-fitted" lines to "badly- fitted" lines. The progress of the algorithm can be assessed 
by monitoring the proportion of well-fitted lines to badly-fitted lines. One would expect 
this ratio to increase as the algorithm progresses. If this ratio decreases significantly it 
may indicate that the line length has been increased too rapidly. By detecting this 
condition, the line length can be safely decreased until the ratio of well-fitted lines to 

15 badly-fitted lines increases. 

After several iterations, the data plot appears as in Figure 7. In this Figure, the whole 
signal appears continuously on the plot with breaks for the phase transitions. The fitted 
line lengths now cross the transitions without undue effect. This is because the data is 
20 centred on the cyclic mean, so the transition points above the trend line have similar 
effects to the points below the trend line. Even if the transitions are asymmetric the gain 
in accuracy resulting from using longer sections more than offsets the loss caused by the 
asymmetry. 

25 Figure 8 shows the final result of the algorithm after the carrier frequency has been 
estimated using the whole signal. The phase versus time plot is now a horizontal line, 
and the carrier frequency f c of the signal is the sum of the estimated carrier frequencies 
f c % that have been removed from the signal. 
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A phase-modulated signal to which the algorithm can be successfully applied must have 
the following properties: 

1. The signals must be phase-modulated continuous signals, where the modulation is in 
the form of quick phase jumps when compared to the rest of the signal; 

2. The approximate size of the phase jumps must be known, and the size must be even 
such that an integer multiple of jumps fits into 2n; 

3. Other than during phase modulation the signals must have nearly constant frequency 
(since we are trying to measure a constant frequency); and 

4. Phase transitions should only make up a small proportion of the total signal - 
approximately 30% is acceptable. 

Some advantages of the method disclosed herein over existing methods used for 
frequency estimation will now be discussed. Firstly, in the method of the invention 
phase transitions do not need to be identified or removed. This is because the cyclic 
mean is unwrapped to the middle of the phase range at each iteration (to n/2 in the case 
of n phase jumps), so approximately the same number of transition points exist in each 
half of the phase range. This reduces the effect of transition points on the fitting of lines 
to the data. More accurate results can be obtained by using a robust least-squares 
method for fitting a straight line to the data. Robust least-squares methods use median 
values as opposed to mean values and are therefore less sensitive to outliers in the data. 

Secondly, when using the method of the invention phase transitions can occur in either 
direction without affecting the results. When a phase jump occurs it is wrapped back 
into the relevant range, so the direction of the jump does not matter. The method will 
work with any integer multiple of phase jump, provided the noise is smaller than the 
size of the phase jump. 

Another advantage of the method is that wrapping is deterministic and robust, in 
contrast with conventional unwrapping where the phase transitions' direction and size 
need to be known or determined somehow during processing. If noise is present it can 
be difficult to calculate these attributes, and so the unwrapped data may have spurious . % 
jumps present that will affect the gradient of the line fitted to it. The proposed wrapping 
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method needs no assumptions to be made about the phase before and after the transition, 
since it does not use this information. It is therefore robust to noise. 

A further advantage of the method is that it is very accurate. This is because all of the 
5 data is used in the final frequency estimate. It has accuracy approaching that produced 
by a conventional line fit to perfectly unwrapped data where the transitions have been 
removed. It is also much more accurate than using an average of many short line fits; 
Conventionally, if transitions cannot be removed only the continuous sections between 
them can be used to estimate the frequency, and these are very short compared to the 
10 total data length. . 

Yet another advantage of the method is that it is iterative and always converges as long 
as a small enough rate of change of line length is chosen. The algorithm iterates to a 
good solution by fitting small lines to the data and then adjusting the data and fitting 

15 longer lines. In this way, the data are gradually all brought into range so that a long, 
continuous line can be fitted to give an accurate estimate of the gradient. If the increase 
in line length / from one iteration to the next is too large the algorithm can diverge. 
However, this can be recognised by examining the accuracy of fit of the lines to the 
data. The accuracy measure can then be used to produce an algorithm with an adaptive 

20 line length that will bring the algorithm back to convergence. 

This method can be used in any application where accurate estimation of signal 
frequency is required, but where the signal is phase modulated such as in 
communications. The nature of the algorithm means that the receiver of the signal need 
25 not know the rate or exact form of the modulation, only the number of phase states. The 
method can also be used to identify Doppler shift in signals such as satellite 
transmissions received in moving vehicles, or to "lock onto" a phase-modulated carrier. 
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